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Abstract 

Conservative algorithms for boundary interfaces of 
overlaid grids are presented. The basic method is ze- 
roth order, and is extended to a higher order method 
using interpolation and subcell decomposition. The 
present method, strictly based on a conservative con- 
straint, is tested with overlaid grids for various applica- 
tions of unsteady and steady supersonic inviscid flows 
with strong shock waves. The algorithm is also applied 
to a multi-level grid adaptation in which the next level 
flner grid is overlaid on the coarse base grid with an 
arbitrary orientation. 

I. Introduction 

Unstructured grid algorithms [1-3] have been pop- 
ular in recent years, because they have flexibility in grid 
topology and also have efficiency in grid adaptation, in 
comparison with structured grid algorithms. However, 
the main drawback of the unstructured grid algorithms 
lies in the limitation of selecting the time integration 
scheme: one has to use an explicit or point- wise implicit 
algorithm. A general motivation of the present study is 
to develop a counterpart of the unstructured grid algo- 
rithms for structured grids. This counterpart must not 
only have the geometric flexibility and grid adaptation 
efficiency of the unstructured grid algorithm, but also 
must have a freedom of choice in a variety of existing 
implicit schemes. 

To be specific, the present methodology use the 
overlaid grid (Chimera technique) [4-5] to improve geo- 
metric flexibility, and the multi-level grid adaptation 
[6,23] for efficient grid adaptation. Multi-level grid 
adaptation here is not the mesh embedding technique 
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in which coarse grid cells are locally refined by embed- 
ding the next-level finer subcells. Instead, a next-level 
structured grid block is generated overlaying on the 
coarse base grid with an arbitrary orientation. This 
approach has the following merits. First, it gives geo- 
metrical independence between grid levels. Second, it 
allows natural alignment and clustering between grids 
and features in the solution. This grid alignment 
with the solution feature (e.g, shock or contact sur- 
face) possesses an extra advantage in conjunction with 
characteristic-based shock capturing shemes such as the 
flux-difference splitting [7] or flux-vector splitting [8] 
schemes. Since these numerical schemes are, in gen- 
eral, based on a one dimensional theory, if a shock or 
contact discontinuity is obliquely aligned with the grid 
cell (typically in a two dimensional case), the disconti- 
nuity cannot be numerically captured in fewer cells as 
when the discontinuity is aligned with the grid. 

Overlaid grid and multi-level grid adaptation will 
introduce blocks of structured grids which can overlay 
each other with an arbitrary orientation. Then each 
block of the structured grid is solved by an implicit 
scheme using an ADI or LU algorithm. This approach, 
which is generally defined as a domain decompostition 
method, is also naturally suited for parallel processing. 

However, the success of the overlaid grid and multi- 
level grid adaptation techniques relies on correct com- 
munication among overlaid grids through boundary in- 
terfaces. Previous studies on this problem can be 
catagorized into two groups: (i) interpolation of the 
neighboring variables [4,5,9], and (ii) characteristics- 
type approach [10-12]. Since neither of them necessarily 
satisfies any form of a conservative constraint, the solu- 
tions on overlaid grids are often mismatched with each 
other. This may result in instablility or oscillations, es- 
pecially when a shock wave passes through boundaries 
of overlaid grids. Nonetheless, characteristics-type ap- 
proaches show improvements in accuracy and conver- 
gence over interpolation methods. 

The main constraint in updating boundary inter- 
face variables is that any valid procedure must ensure 
conservative properties. For the generally overlaid grid 
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system, Berger [13] attempted to derive a conservative, 
finite difference flux scheme. However, the stencil of 
fluxes becomes irregular and consequently raises diffi- 
culties in its extension, not only to three dimensions 
but also to implicit implementations (e.g, extra Jaco- 
bians emerge with consequent restructuring of the left 
hand side matrix). In the early work, Rai [14] applied a 
conservative procedure to patched grid systems where 
the interface between two patch boundaries becomes a 
line for two dimensions. In this case global conservation 
can be implemented as follows: balance of the temporal 
surface fluxes and continuity requirment of the depen- 
dent variables at the patch interface. This approach has 
been extended to three-dimensional patched grids [17- 
19] where the interface between grid patches becomes a 
plane. The basic idea is to obtain the unknown surface 
fluxes by interpolation with the area-weighted coeffi- 
cients, enforcing a global conservative constraint. 

In the present study with overlaid grids, we ap- 
ply a conservative constraint to an arbitrary domain in 
the region of overlap such that the conserved quanti- 
ties such as the mass, momentum, and total energy are 
preserved in the domain at a particular time. Similar 
ideas have been used in the Lagrangian hydrodynamics 
[15,16] for the rezoning process (transfer of solutions 
from an old grid to a new grid). Under the conser- 
vative constraint, we interpolate the unknown depen- 
dent variables Q in the domain by area-weighting. We 
also extend the basic method (zeroth order) to a higher 
order algorithm by reconstructing the given piecewise 
constant solution, cell-by-cell, under the restriction of 
conservative constraint (volumetric conserved quanti- 
ties in a cell). Application of the higher order method 
becomes important, especially when the data are trans- 
fered from coarse cells to fine cells. 

In the following sections, we describe the governing 
equations, basic solver algorithms, and conservative al- 
gorithms for general overlaid grids. A number of appli- 
cations are presented for the unsteady and steady in vis- 
cid problems with overlaid grids. Finally, this conser- 
vative procedure is applied to a multi-level grid adap- 
tation. 


II. Governing Equations and Solver Algorithms 


The integral form of the time-dependent compress- 
ible Euler equations can be written in two dimensions 


TtJI A QdA + I s ™ ds = ° 


( 1 ) 


Q = (p, pu, pv y e) T is a vector of conservative variables 
and the vector F denotes the flux of mass, momen- 
tum, and energy across the surface of the control vol- 


ume. Here e represents the total energy per unit vol- 
ume. Based on the finite volume method, equation (1) 
is semi-discretized by assuming that the cell-centered 
conserved variables are constant within a cell, and that 
the flux integral at cell surfaces is also approximated by 
an average value of the numerical flux and the surface 
length. 

We follow the MUSCL approach and Van Leer’s 
flux-vector splitting procedures in evaluating the nu- 
merical flux at the surfaces. To be specific, we first use 
upwind-biased interpolation of a certain set of variables 
to obtain the values at the cell surface. With the in- 
terpolated values, the split fluxes at the cell surface are 
then evaluated according to flux-vector splitting pro- 
cedures. The scheme is fully upwinded, second order 
accurate or upwind-biased, third order accurate for the 
convective and pressure terms. More details on the 
scheme can be found in reference [8]. Here we would 
like to briefly remark on the convergence of the scheme. 
We found that convergence is sensitive to the choice of 
the interpolating variables in the MUSCL approach and 
to the choice of limiters as well. The best convergence 
is achieved with primitive variables and the diffusive 
min-mod limiter. In the case of conserved variables 
for interpolation and the compressive min-mod limiter, 
convergence is locked in a limit cycle. 

In the present study, the fourth order Runge-Kutta 
time-stepping scheme is used for the unsteady calcula- 
tions. For the steady state problems, the first order 
Euler implicit method is used with an ADI algorithm. 


III. Conservative Algorithms for Overlaid Grids 

Suppose we have two domains which are overlaid 
in a general manner (see figure 1). Let Q 1 and Q 2 
represent sets of conservative variables for the domains 
1 and 2, respectively. A particular section in the region 
of overlap is denoted by a domain O. 

At a particular time t, a conservative constraint 
over the domain Q is represented by 

f Ql(x,t)dQ = f Q2(x,t)dQ = C(t) (2) 

Jci Jft 

where Q is a function of time t and space coordinates 
x. For the rest of this section, let a distribution of Q 2 
over the domain Q be known. Then a set of conserved 
quantities C(t) (e.g, mass, momentums, and total en- 
ergy) over the domain Q is known at time t, according 
to the equation (2). In a discrete sense, two domains 
are represented by computational grids and the piece- 
wuse constant, discrete solution of Q 2 over the domain 
Q is given. 
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The basic idea is to obtain the unknown Ql in 
the domain Cl from the given Q 2 by interpolation with 
the area-weighted coefficients, enforcing the conserva- 
tive constraint, equation (2). There are several impor- 
tant facts to keep in mind. First, equation (2) is only 
a constraint condition over the domain Cl. Therefore 
we would like to apply the equation (2), cell-by-cell, so 
that the conservative constraint is satisfied automati- 
cally in a global sense. Second, the constraint condi- 
tion concerns only the conserved quantities C(t) over 
the domain Cl (e.g, a Q2 cell), not the variation of the 
conservative variable Q 2(x, t) within the domain. In 
other words, the conservative constraint is satisfied, as 
long as the conserved quantities are preserved over the 
domain Q. This leads to a freedom in reconstructing 
the variation of Q 2 within the cell for the purpose of 
achieving a higher order redistribution to Ql. Here we 
designate redistribution as the process of determining 
the distribution of the unknown Ql over the domain 
O, under the conservative constraint. Third, the ac- 
curacy of conserved quantities C(t) obtained for the 
domain Cl was predetermined by the truncation error 
in the solution of the grid 2. This means that the ac- 
curacy of the transferred conserved quantities on grid 

1 cannot be increased by increasing the order of the 
redistibution procedure. Consequently, the best pos- 
sible solution which can be constructed on grid 1 by 
this conservative constraint is the solution of the grid 

2 itself. 


Zeroth Order Redistribution 


Suppose the domain Cl corresponds to a single cell 
of coarse grid 1 and to composite pieces of cells of fine 
grid 2. A schematic is shown in figure 2. Given a 
piecewise constant distribution of Q2 from the solution 
on grid 2, apply the conservative constraint over the 
domain O. Then equation (2) yields, in a discrete form, 




Q2i t i 


i = l L 


)/,r 


(3) 


where (jf, k) and (l, m) represent cell indices for grid 1 
and 2, respectively. Also index i indicates each (/, m) 
cell overlapped by the (j, k) cell, and imax is the total 
number of them. Aj t * is the area of the (j, k) cell, and 
(i OT )/ |m is the portion of the (/, m) cell area overlapped 
by the (j, k) cell. This situation corresponds to the case 
of transfering the data from the fine cells to a coarse 
cell. The conservative constraint itself is sufficient to 
uniquely determine the Ql variables. We call this pro- 
cedure the zeroth order redistribution, since we use the 
given piecewise constant values of Q 2. 

Conversely, suppose the domain Cl corresponds to 
cells of fine grid 1 and to a single cell of coarse grid 2 


(see figure 3). This case corresponds the data transfer 
from a coarse cell to fine cells. Then the conservative 
constraint yields again, in a discrete form, the following: 

imax 

yi [Q^j,h(^ov)j,k] = Q2j >fn A/ jm i E (j, k) (4) 

izzl 

where index i indicates each (jr, k) cell overlapped by 
the (/, m) cell, and (A^j^ represents the portion of 
the (j, k) cell area overlapped by the (l, m) cell. In 
this situation, the conservative constraint is only a nec- 
essary condition, but is not sufficient to uniquely de- 
termine the Ql variables (one equation with imax un- 
knowns). The zeroth order redistribution will set Ql 
equal to Q 2 over the domain Cl. Then equation (4) is 
automatically satisfied and the conservative constraint 
still holds. However, it yields an undesirable property 
such that the redistributed finer cells will represent a 
continous distribution of the coarse cells as a stairstep 
approximation (i.e, degrading smothness). This may 
have effects on stability and consequently on the con- 
vergence. 

Higher Order Redistribution 

In order to circumvent this problem, a higher order 
redistribution procedure is proposed here. Since equa- 
tion (4) is an underdetermined system, we have free- 
dom to assign a proper value to each overlapping (j, k) 
cell, while a conservative constraint is maintained. The 
idea is based on a reconstruction of the given piecewise 
constant solution over a Q2 cell with the conservative 
constraint. 

It would be nice if we could represent the piecewise 
variation in a single, analytic form, with the conserva- 
tive constraint still satisfied in a Q 2 cell (i.e., preserving 
the volumetric conserved quantities). For example, a 
linear approximation for the reconstructed Q 2 within a 
cell is expressed as an equation of a plane with a single 
free parameter which is determined by the conservative 
constraint. However, the only fact that matters is how 
many cells of grid 1 are occupying the grid 2 cell (i.e, 
area ratio between two grid cells), because, after all, the 
occupying grid 1 cells will be represented as a piecewise 
constant approximation. Therefore we propose to rep- 
resent a variation within a Q2 cell with a finite number 
of piecewise constant approximations. 

If we have a coarse (l, m) cell of grid 2, then de- 
compose it into multi-subcells denoted by index is. For 
example, decompose the (l, m) cell into four subcells by 
connecting the midpoints of four surfaces of (/, m) cell. 
A schematic is shown in figure 4. These four subcells 
are not necessarily to be matched with the fine cells of 
grid 1 in figure 3, but it is important to keep the area 
ratio between grid 1 cells and the decomposed subcells 
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close to unity. The best strategic decomposition is to 
mimic the orientation and size of the cells of grid 1. 

Now obtain an interpolated value Qi s for each sub- 
cell (/, m, is). For example, Qu - 4 is a subcell-centered 
value obtained by triangular interpolation using its sur- 
rounding coarse cells, t.e, <?2j >m , C?2j_i )Tn , and 
(see figure 4). Then we define a weighting coefficient 
Cu for each (/, m, is) subcell as follows 


C u = 




£i/f=l Q id Aid 


is £ (1, m) 


(5) 


where Qis — Q**(Q2i >m ), A% 9 is a subcell area, and id 
represents a dummy index for subcell index is. The 
weighting coefficient is designed to be proportional to 
the product of the interpolated value and its area so 
that it relates to the conserved quantities, while it is 
scaled so that summation of all Cu becomes unity. 

Multiplying the volumetric conserved quantities 
Q2i tTn Ai >m by the weighting coefficient Ci Si finally re- 
constructs the piecewise constant Q 2 over the (/, m) 
cell into four piecewise constant Qi $ over the (/, m, zs) 


subcells where Q is is expressed as 


l^id=l QidAid 


( 6 ) 


Then the reconstructed subcell value Qi s automatically 
satisfies the conservative constraint over the domain Q 
(i.e, the (1, m) cell) as 


4 

\Qi$Ai$\ = Q2l,mAl } rn (7) 

is = 1 

This concludes the process of reconstructing Q 2 over 
the domain H for the four subcell case. 

We now generalize the above procedure to the case 
where the total number of subcells within the (/, m) cell 
is ismax . Returning to the case described in figure 3, 
equate the left hand sides of equations (4) and (7) to 
obtain 


imax ismax 

E ~ E- \Q%»Aig\ (8) 

t—l iszzl 

If ismax > imax , then the situation becomes such that 
the conservative constraint can uniquely determine the 
Q 1 variables: 


< ismax 


Qij.k = Y2 


iQ; (^)* 


is=l 


Aj t k 


is € (/, m) (9) 


Finally, the conservative constraint, equation (2), 
can be written in a general form for the higher order 
redistribution as 



where i indicates each overlapped (/, m) cell and is in- 
dicates each overlapped subcell (1, m, is) within (I, m) 
cell. The order of accuracy in this redistribution is de- 
pendent upon the selected order in the interpolation 
procedure and the total number of subcells, ismax , 
within an (/, m) cell as well. However, if the upper 
limit of is equals to ismax in equation (10), then the 
equation (10) automatically becomes equation (3). 

The reconstruction process described here can be 
viewed as the procedure of raising the order of accuracy 
of the initial-value interpolation in a “projection stage” 
[20] of a Godunov scheme. The higher order redistribu- 
tion may create a non-physical solution (overshoot or 
undershoot) near discontinuities. Therefore a limiter 
is recommended for future application in the interpo- 
lation process. For all the test cases presented in this 
paper, the limiter was not necessary. 

Algorithm for Computing Area of Overlap 

The conservative algorithms described in the previ- 
ous sections can be summarized as a rule of area weight- 
ing. The remaining question is how to evaluate the area 
of overlap between two arbitrary shaped quadrilateral 
meshes. The algorithm should consist of logically sim- 
ple procedures that are extendable to the three dimen- 
sional case (evaluating a volume of overlap). 

The algorithm consists of two processes called “en- 
closure” and “crossing”. “Enclosure” determines if a 
point is within a cell or not, and “crossing” calcu- 
lates the crossing point between either a line with fi- 
nite length or one with infinite length. Using these 
two processes, one can find the vertices of the area of 
overlap, which always form a convex polygon. The pro- 
cedure consists of six possible combinations in forming 
the vertices of a polygon from the vertices and crossed 
points between two overlapping quadrilateral meshes. 
Finally the vertices of the polygon are arranged in a 
connecting order. 

The area of the convex polygon is evaluated by 

1 n r 

A °v - o I - x iyi+i ) I ( n ) 

*=1 

where np represents the total number of polygon sur- 
faces. Here (z;, t/;) denotes a vertex point of a polygon. 
The most computationally intensive part of this algo- 
rithm lies in determining which cells of two grids share 


where (Awji* represents the portion of (/, m } is) sub- 
cells overlapped by the (j, k) cell. 
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a common region. A knowledge-based tree search al- 
gorithm was used to improve the efficiency. This is a 
pre-processing procedure with a cost of memory over- 
head for the geometric information. 

IV. Results and Discussion 

The conservative algorithm for overlaid grids, de- 
scribed in the previous section, was applied to vari- 
ous unsteady and steady problems of supersonic in- 
viscid flows. In order to demonstrate the conserva- 
tive property of the present method, all the problems 
possess strong moving or stationary shock waves which 
cross the boundary interfaces of the overlaid grids. Fi- 
nally, this procedure was combined with multi-level grid 
adaptation to display the versatility of the algorithm. 

Before discussing the results, some details of im- 
plementing the conservative interface treatments will 
be briefly described. Since the fully upwind second or- 
der or upwind-biased third order accurate scheme has 
a five point stencil, the scheme reduces to first order 
accuracy at the two cells next to a boundary. There- 
fore two cells at the boundary of grid 1 are intentionally 
placed away from the other two cells at the boundary 
of grid 2 (shown in figure 5). There are two reasons 
for this. First, the third cells from the boundaries of 
grid 1 and 2 can use the five point scheme to have at 
least second order accuracy. Second, the two boundary 
cells of grid 1 and 2 can be updated by the conser- 
vative algorithms using at least second order accurate 
solutions from each grid. This boundary update is ex- 
plicitly conducted so that it can be easily applied to 
explicit or implicit schemes. 

One last remark regards the synchronization of the 
solutions on both grids. The conservative constraint, 
equation (2), is valid at a particular time on both grids. 
Therefore, if two different time steps are used for each 
grid, one must ensure that the conservative boundary 
update is conducted at the same elapsed time, espe- 
cially for the time accurate calculations. One way to 
synchronize is to interpolate the solution in time. For 
example, advance a solution on one grid with large time 
step. Then, interpolate the solution in time to a num- 
ber of time levels necessary for the other grid solution 
to “catch up” to the advanced time. Finally apply the 
conservative algorithm at each sub- time level, using the 
timewise interpolated solution. The accuracy and sta- 
bility of this procedure is not yet tested. For all the re- 
sults shown in the present study, synchronization was 
achieved by advancing the solutions with the smaller 
time step between the two grids. 

Time-dependent Calculations 


Case 1 

The conservative algorithm is tested for a time- 
dependent, shock tube problem. Here we intention- 
ally solve this one-dimensional problem in two dimen- 
sions. The computational domain is composed of over- 
laid grids: 45x9 and 55x11 meshes for grid 1 and 2, 
respectively (shown in figure 6 (a)). The lines between 
two overlaid grids are mismatched both horizontally 
and vertically. Pressure and density ratios of 10 are 
set at x = 2 as an initial condition. 

Figure 6 (b) displays the density contours at 
t=0.61, which are plotted on top of each other. The 
zeroth order redistribution matches the solutions well 
in the overlapped region, since the grid spacings of the 
two grids are comparable. A comparison of the so- 
lutions with the exact solution is shown in figure 6 
(c). The solid line represents the exact solution (den- 
sity distribution), and the open circles and the dashed 
line correspond to the solutions from grid 1 and 2, re- 
spectively. Results show that the data communication 
through boundary interfaces between two overlaid grids 
was conducted conservatively and time-accurately. 

Case 2 

The same problem as case 1 is tested, but the 
grid spacing between the two grids is doubled in the 
x-direction: 41x9 and 81x11 meshes for grid 1 and 2, 
respectively (shown in figure 7 (a)). Results of the ze- 
roth and higher order redistributions are compared in 
figure 7 (b) and (c). Figure 7 (b) (the zeroth order 
redistribution) shows the mismatched density contours 
from the two solutions in the overlapped region. Ac- 
tually, the coarse grid solution is accurately obtained, 
because the data are transfered from fine cells to coarse 
cells. The difficulty is encoutered for the fine grid in 
which the situation is reversed. The stairstep contours 
can be observed in the figure 7 (b). The result of the 
higher order redistribution, shown in figure 7 (c), no- 
ticeably improves the solution of the fine grid in the 
overlapped region. For the results from the higher or- 
der redistribution presented in this paper, the coarse 
cells are decomposed into four subcells, according to 
midpoints of the coarse cell surfaces. As mentioned 
earlier, a more optimal subcell decomposition will lead 
to a better match of the solutions between two overlaid 
grids. 

Case 3 

A more realistic two dimensional unsteady inter- 
actions of a moving shock of M t = 7.19 in a 20° bent 
channel were simulated on two overlaid grids, using the 
conservative algorithm with the higher order redistri- 
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bution. Figure 8 (a) shows the two overlaid grids with 
101x27 and 91x31 meshes. A moving normal shock of 
M, = 7.19 is set at the inlet of the 20° bent channel as 
an initial condition. 

This bent channel consists of convex and concave 
corners where expansion and compression of the in- 
duced supersonic flows (M = 1.7721) occur, respec- 
tively. In the lower compression region, a complex- 
Mach reflection, which is defined as a single Mach re- 
flection with a kink near the triple point, is developing. 
In the upper expansion region, a Prandtl-Meyer fan is 
centered at the convex corner, and a normal shock wave 
is embedded between the fan and the incident shock 
wave to adjust the flow speed to the condition behind 
the diffracted part of the incident shock wave. Sequen- 
tial plots of density contours (figures 8 (b) - (f)) display 
developments of the interactions among shock waves, 
expansion waves, and confined geometries. The solu- 
tions of two grids are well matched in the overlapped 
region. This result shows the capability of the present 
conservative algorithm for strong multiple shock waves 
crossing the boundaries of the overlaid grids. 

Steady-state Calculations 
Case 1 

As a steady state application, the conservative al- 
gorithm was tested for a two dimensional oblique shock 
reflection on a flat plate. The free stream Mach num- 
ber Mqo is 2.9 and the incident oblique shock angle is 
29°. A uniform mesh of 61x17 is used for grid 1, while 
grid 2 has a uniform mesh of 71x15 (figure 9 (a)). Fig- 
ure 9 (b) shows a close-up view near the overlapped 
region. Mesh cells of each grid are mismatched in both 
x and y directions. Since the grid spacings between 
the two overlaid grids are comparable, the zeroth order 
redistribution was used. Figure 9 (c) displays pressure 
contours of each grid plotted on top of each other. The 
solutions from the two grids are well matched in the 
overlapped region. Two parallel solid lines indicate the 
boundaries of the overlapped region. The performance 
of the conservative algorithm for steady state problems 
is expected, since we iterate the unsteady solver to the 
steady state in a time marching manner. 

Case 2 

In this case, the incident oblique shock is simu- 
lated in the domain, 0 < x < 1.6 and 0 < y < 1. with 
the same physical flow conditions as the previous case. 
This time, the Rankine-Hougoniot conditions are ap- 
plied at (0, 0). A uniform mesh of 41x17 is used for 
grid 1, while grid 2 has a uniform mesh of 81x15. The 
close-up view of the overlapped region is shown in fig- 


ure 10 (a). The grid spacing ratio between the two 
grids is 2 : 1 in the x direction. Figure 10 (b) shows the 
pressure contours of the zeroth order redistribution. As 
discussed in the previous comparison for the unsteady 
case, stairstep contours of the fine grid can be observed 
in the lower portion of the overlapped region. Figure 

10 (c) displays the pressure contours improved by the 
higher order redistribution. 

Case 3 

Supersonic flow with a free stream Mach number of 
Mqo =6.0 over a blunt cylinder is tested with two over- 
laid grids. Grid 1 of 71x30 is placed next to the body, 
and grid 2 of 81x35 is overlaid on the grid 1 (shown in 
figure 11 (a)). Figure 11 (b) displays the details of grid 
lines in an upper portion of the overlapped region. The 
zeroth order redistribution is used in this case. Figure 

11 (c) shows Mach number contours plotted on top of 
each other. The two solid lines indicate the boundaries 
of the overlapped region. Mach number contours of ex- 
panding flow from low subsonic to moderate supersonic 
are well matched in the overlapped region between the 
two grids. The surface pressure distribution is also com- 
pared with experiment in figure 11 (d), and agrees well 
with it. 

Case 4 

Edney’s type III shock-on-shock interaction on a 
blunt cowl lip was tested with two overlaid grids. Su- 
personic flow with a free stream Mach number of M^ = 
6.0 over a blunt cowl lip is intersected at 9 — 14.4° by 
an oblique shock inclined to the free stream at an an- 
gle of 22.67°. More details on the flow conditions can 
be found in [21]. The grid 1 of 89x39 is placed next 
to the body, and grid 2 of 99x49 is overlaid on grid 1 
(figure 12 (a)). Figure 12 (b) shows a close-up view of 
the upper portion of the overlapped region. For this 
case, the zeroth order redistribution is also used. Fig- 
ure 12 (c) shows Mach number contours plotted on top 
of each other. The overlapped region encloses interact- 
ing shocks, a triple point, shear layer, and expanding 
flows. Two solid lines indicate the boundaries of the 
two overlaid grids. Mach number contours of the two 
overlaid grids are well matched. The surface pressure 
distribution is also compared with experiment in figure 

12 (d), and agrees relatively well, except near the shear 
layer impingement region. Better agreement would be 
obtained not only by solving the Navier-Stokes equa- 
tions with the proper turbulence model for this turbu- 
lent shear layer [22], but also by refining the meshes of 
the two grids in the region near the impinging shear 
layer. 



Multi-level Grid Adaptation 

The present conservative algorithm has a special 
application to multi-level grid adaptation for achieving 
high accuracy. A major difference from the work of 
Berger [23] is the use of next level finer grids which can 
be overlaid on the base grid with an arbitrary orienta- 
tion. 

Briefly describing the adaptation procedures, first 
obtain a solution converged to a certain level (e.g, tru- 
cation level) on a relatively coarse grid. Then, based on 
an error estimator or scanning of the flow with a proper 
detection criterion, generate a next level grid which is 
clustered and aligned with the solution features such 
as a shock wave, contact surface, boundary layer, or 
shear layer. On this grid, the coarse grid solution is in- 
terpolated by using the conservative algorithm. Then 
solutions of the base and next level grids will be ad- 
vanced by communication of data through boundary 
interfaces. 

As a demonstration of the above approach, case 
2 in the steady state calculations was selected for grid 
adaptation. Figure 13 (a) shows a single uniform grid of 
41x41 for a domain (0 < x < 2.5 and 0 < y < 2). Pres- 
sure contours are shown in figure 13 (b). Figure 13 (c) 
displays a pressure distribution at y = 0.5. The oblique 
shock is numerically captured within appoximately 7 to 
8 cells. This is a typical result of Van Leer’s flux- vector 
splitting MUSCL scheme (upwind-biased third order 
accurate) when a shock wave lies obliquely to the grid 
lines, since the scheme is based on the one-dimensional 
theory. 

A next level grid of 33x15 mesh is generated by 
scanning the maximum pressure gradient points from 
the base grid solution. Figure 14 (a) shows the grid 
which is clustered and aligned with the shock. Pres- 
sure contours from the refined solution obtained on the 
next level adapted grid is shown in figure 14 (b). Figure 
14 (c) displays the pressure distributions at three dif- 
ferent locations. As predicted [8] for a one dimensional 
case, the shock is numerically captured within two cells 
for all three locations. This grid alignment is one way 
of achieving the full quality of upwind TVD schemes, 
whereas the rotationally biased upwind scheme [24] has 
shown a similar performance. 

V. Conclusions 

Conservative algorithms for boundary interfaces of 
overlaid grids have been presented. The test results 
show the potential of the present method in various 
applications for unsteady and steady supersonic flows 
with strong shock waves. The proposed higher order 


redistribution presents a distinct improvement over the 
zeroth order redistribution for cases where the data are 
transfered from coarse cells to fine cells with the ratio 
of grid spacing between two overlaid grids no less than 
2:1. The present conservative algorithm has a di- 
rect application to the domain decomposition method 
in which overlaid grids and multi-level grid adaptation 
are currently considered. It is believed that implemen- 
tation of parallel processing for this multi-block domain 
strategy is an immediate next step. 
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Fig. 1 Schematic of two overlaid grids 


Fig. 4 Schematic of subcell decomposition and triangular interpolatioi 
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Fig. 2 Schematic of a uniquely determined system (fine to coarse) 


Fig. 5 Schematic of boundary interfaces layout 
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Fig. 3 Schematic of an underdetermined system (coarse to fine) 
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Fig. 6 (a) Shock tube: grid 1 (45x9), grid 2 (55x11) 
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Fig. 7 (a) Shock tube: grid 1 (41x9), grid 2 (81x11) 
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Fig. 6 (b) Density contours (zeroth order redistribution) 
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Fig. 7 (b) Density contours (zeroth order redistribution) 



Fig. 7 (c) Density contours (higher order redistribution) 


Fig. 6 (c) Density distribution ( — :exact, o:grid 1, :grid 2) 




Fig. 8 (a) 20° bent channel: grid 1 (101x27), grid 2 (91x31), (b)-(f) Density contours of moving shock 
interactions at five sequential stages (higher order redistribution). 
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Fig. 9 (a) Oblique shock reflection: 

grid 1 (61x17), grid 2 (71x15) 


Fig. 10 (a) Oblique shock: Close-up of region of overlap, 
grid 1 (41x17), grid 2 (81x15) 
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Fig. 9 (b) Close-up of region of overlap 


Fig. 10 (b) Pressure contours (zeroth order redistribution) 


Fig. 9 (c) Pressure contours (zeroth order redistribution) Fig. 10 (c) Pressure contours (higher order redistribution) 
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